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ABSTRACT 

We calibrate here cosmological radiative transfer simulations with 
ATON/RAMSES with a range of measurements of the Lya opacity from QSO 
absorption spectra. We find the Lya opacity to be very sensitive to the exact timing 
of hydrogen reionisation. Models reproducing the measured evolution of the mean 
photoionisation rate and average mean free path reach overlap at z ~ 7 and predict an 
accelerated evolution of the Lya opacity at z > 6 consistent with the rapidly evolving 
luminosity function of Lya emitters in this redshift range. Similar to ” optically thin” 
simulations our full radiative transfer simulations fail, however, to reproduce the 
high-opacity tail of the Lya opacity PDF at 2 > 5. We argue that this is due to 
spatial UV fluctuations in the post-overlap phase of reionisation on substantially 
larger scales than predicted by our source model, where the ionising emissivity is 
dominated by large numbers of sub-L* galaxies. We further argue that this suggests a 
significant contribution to the ionising UV background by much rarer bright sources 
at high redshift. 

Key words: Cosmology: theory - Methods: numerical - diffuse radiation - IGM: 
structure - Galaxy: evolution - quasars: general 


1 INTRODUCTION 

The Lya forest is a valuable probe of the underlying mat¬ 
ter density field, as well as of the temperature and ionisa¬ 
tion state of the Inter-Galactic-Medium (IGM) at high red¬ 
shift (see |Rauch|p~998 Meiksin|[2009 for reviews). There is 
strong evidence that, despite a significantly increased effec¬ 
tive Lya optical depth, even the highest redshift QSOs at 
z ~ 6 still probe the post-reionisation highly-ionised IGM 
( Bolton V Haehnelt|2007 cf. Lidz et al.|2 006 and Mesinger 


2010), apart from perhaps the notable z — 7.085 QSO ULAS 


J1120+0641 which appears to show signs of a red damping 
wing suggesting a volume fraction of neutral hydrogen of 

2006} 


10% or more (Fan et al. 2006[ Mortlock et al.||2011 and 
|Bolton et al. 2011), but see Bosman V Becker (2015) (in 
prep) for a recent reassessment of the significance of the 
weak Lya emission of ULAS J1120+0641. 


This is in good agreement with the constraints on the 
reionisation history by the CMB which also suggests that 
most of hydrogen in the Universe has been reionised earlier 
than z ~ 6, even though the most recent Planck (Planck Col¬ 


laboration et al. 2015) measurements suggest a somewhat 
later (end of) reionisation than earlier CMB measurements. 

The accelerating evolution of the Lya optical depth at 
z > 5 ( Fan et aL| |2006) and the apparent rapid decrease of 


the observed Lya emission of high-redshift galaxies (Bolton 
V Haehnelt 2013 and Dijkstra et al.|2014 ) also suggest that 


at z ~ 6 we are witnessing the final stages of the reionisation 


Pentericci 


et al. 

2011 

ITreu et al.| 2013, |Caruana et al. 2014[ |Faisst 

et al. 

2014 

[Pentericci et al. 2014, Konno et al.| 2014 and 


by a recent accurate measurement of Lya optical depth PDF 


at 4 < z < 6 by Becker et al. (2015) based on a large sample 


of high-quality high-redshift QSO absorption spectra. 

Simulation of the hydrogen Lya forest at z < 5 have 
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mostly been performed assuming that hydrogen is already 
highly ionised, that the Universe can be assumed to be opti¬ 
cally thin for hydrogen ionising radiation and that radiative 
transfer effects can be neglected. The UV background in 
these simulations is assumed to be spatially homogeneous 
as in widely used spectral synthesis models of the UV back¬ 


ground (Haardt & Madau 1996 Haardt & Madau 2001 
Haardt L Madau|2012). 


Such simulations have been used widely to infer the 
hydrogen photoionisation rate by rescaling the optical depth 


et al. 1997} 

Theuns et al. 1998, Bolton et al.|2005[ Bolton 

& Haehnelt 

2007, Faucher-Giguere et al. 2008a|, Faucher- 


Bolton ]2013[ ). When comparing such simulations to their 


measurements of the Lya optical depth PDF Bec ker et al.| 


(2015) reported rapidly increasing deviations for the high 


opacity tail of the PDF at z > 5 and interpreted that being 
due to spatial fluctuations in the photoionisation rate in the 
post-overlap phase of reionisation as expected in models for 
inhomogeneous reionisation (Miralda-Escude et al. 2000). 


More detailed modelling of this deary requires full cos¬ 
mological radiative transfer simulations. Most cosmological 
radiative transfer simulations have concentrated on the large 


scale topology of reionisation (JGnedin & Abel 

2001, 

Ra- 

zoumov et al. 2002[ Iliev et al. 2006, Lidz et al. 

2007 

Trac 

& Cen 2007[ Aubert & Teyssier 2008, Finlator et al. ! 

2009, 

Petkova & Springel|2009, Aubert & Teyssier 2010, Chardin 

et al. 2012) in rather large simulation boxes. Fully under- 


standing the growth of ionised regions requires, however, 
not only an accurate modelling of the large scale distribution 
of the sources of ionising radiation but also of the sinks of 
ionising radiation as is being increasingly real ised (jMiralda- 


Escude et al.| 2000[ |Furlanetto &; Oh| [2QQ5| |Furlanetto V 
Mesinger|2009 ]Choudhury et al.| 2009[|Mc Quinn et al.|201l[ 


Kaurov V Gnedin] |2013[ |Gnedin| 2014[ Gnedin V Kaurov| 

2014[|Choudhury et al.|2014|and|Mesinger et al.|2015|). 

This makes accurate cosmological radiative transfer 
simulations numerically extremely demanding in terms of 
dynamic range and hybrid techniques to treat small and 
large scale separately are being develope d (|Mesinger V 


Furlanetto|2007| |Choudhury et al.|2014| and |Mesinger et al. 
2015) . 


Accurate modelling of the sinks of ionising radiation is 
particularly important for modelling of Ly<a forest data and 


this is the main focus of the paper here (see also Petkova & 
Springel 2009). For this, we use here a suite of rather high- 


resolution hydro-dynamical RAMSES simulations (Teyssier 


2002) post-processed with the GPU-based radiative trans¬ 


present prospects concerning further work in section [4] Our 
simulations assume the following cosmological parameter (as 
derived from the 2013 Planck temperature power spectrum 
data alone, Planck Collaboration et al. 2 014) : Q m = 0.3175, 
Qa = 0.6825, Q b = 0.048, h = 0.6711, a 8 = 0.8344, and 
n s = 0.9624. 


2 NUMERICAL SIMULATIONS 

We briefly describe here the codes used for both modelling 
the hydrodynamics of the gas and the radiative transfer cal¬ 
culation performed as a post-processing step on the hydro- 
dynamic simulation outputs. 


2.1 Cosmological hydrosimulations with RAMSES 

We have performed cosmological hydro-simulations with the 


adaptive mesh refinement code RAMSES (Teyssier 2002). 


In RAMSES the gas dynamics is solved with a 2nd order 
Godunov scheme combined with a HLLC Riemann solver. 
The collisionless dynamics of the dark matter (DM) is rep¬ 
resented by DM particles and the gravitational potential is 
calculated with a multi-resolution multi-grid solver. Initial 
conditions were produced with MUSIC ( Hahn V Abel|2013 ). 

Star formation is included following the prescription of 


Rasera & Teyssier (2006). The star formation recipe assumes 
that above a baryon over-density of 6 ~ 50 gas transforms 
into stars of constant mass with efficiency e = 0.01 per free 
fall time. We do not include stellar or AGN feedback and 


metal enrichment has also not been included (see Bauer et al. 
|2015| for a recent cosmological radiative transfer simulation 
based on the Illustris simulation with AGN and stellar feed¬ 
back tuned to form a realistic galaxy population). 

The usual primordial cooling processes for hydrogen 
and helium are taken into account: collisional ionisation 
cooling, recombination cooling, dielectronic recombination 
cooling, collisional excitation cooling, bremsstrahlung and 
inverse Compton cooling (see T heuns et al.| 1998). Simula¬ 
tion outputs are generated every 40 Myrs from z ~ 30 to 
z ~ 2 which results in 80 outputs for each of the simulations 
described further in section [2~4l 


2.2 Radiative transfer in post-processing with 
ATON 

The radiative transfer calculations have been performed as 
a post-processing step with the ATON code. ATON is de¬ 


fer code ATON (Aubert & Teyssier 2008). Note that the scribed in Aubert & Teyssier (2008) and has been mainly 


box size of our simulation are smaller than those of many 
other reionisation simulations in order to properly resolve re¬ 
gions optically thick to ionising radiation in the post-overlap 
phase. 

The paper is organized as follows: In section [2] we 
present the main features of the simulations performed in 
the context of this work. In section [3] we present our re¬ 
sults first for the evolution of global quantities in the sim¬ 
ulation and second for the evolution of Lyman-alpha for¬ 
est statistics in our different model. We discuss our results 
and make predictions for the evolution of the ionising UV 
background in section |T2| and we give our conclusions and 


used to model hydrogen reionisation (Aubert V Teyssier 
Chardin et al.||2012 and Chardin et al.||2014 ) and 


2010 


particular the reionisation of our Local Goup of galaxies 
( Ocvirk et al.|20l3 and Ocvirk et al.|2014 ). ATON employs 
a moment-based description of the radiative transfer equa¬ 
tion, based on the Ml approximation that provides a simple 
and local closure relation between radiative pressure and 
radiative energy density. The code takes advantage of GPU 
acceleration to solve the equations explicitly in conservative 
form while satisfying a very strict Courant condition, i.e. 
At < c Ax/3, where At is the time step, Ax the cell size 
and c the speed of light. 
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Note that the radiative transfer in ATON uses the ac¬ 
tual speed of light c and not a reduced speed of light approx¬ 
imation (see |Gnedin fc Abel 2001 and Bauer et al. 2015). 
Note further that we have not used here ATON’s ability to 
follow photons with a range of frequencies. We have used 
instead the computationally cheaper options of tracing only 
monochromatic ionising photons with energy of 20.27 eV, 
approximately the mean photon energy of a 50000 K black- 
body spectrum (see also Baek et al. 2009). 

The radiative transfer is performed on a grid with reso¬ 
lution equal to that of the coarse base grid of our RAMSES 
simulations. Note that ATON is not able to follow the AMR 
structure of RAMSES and that the RAMSES simulations 
used do not include any level of mesh refinement. Note fur¬ 
ther that recently a version of this radiative transfer scheme 
has been directly implemented in RAMSES (Rosdahl et al. 


2013) that couples self-consistently radiative forces to the 


the dynamics of the gas. This version is, however, compu¬ 
tationally much more expensive and we do not use it here. 
During post-processing the density and temperature fields 
from the RAMSES simulations are updated within the ra¬ 
diative transfer calculation every 40 Myrs corresponding to 
the frequency of hydrodynamic simulation outputs obtained 
with RAMSES. 


2.3 Modelling of ionising sources 

2.3.1 Optically thin simulations with the HM2012 UV 
background model 

Extensive previous modelling of Lya forest data by some of 
the authors was based on hydrodynamical simulations with 
P-GADGET3 (an improved version of GADGET2 last de¬ 
scribed in Springel 2005) neglecting radiative transfer. To 


make contact to this work we have first compared “optically 
thin” RAMSES simulations with a homogeneous UV back¬ 
ground with comparable P-GADGET3 simulations to make 
sure that we obtain consistent results. 

For this purpose, we have implemented the latest ver¬ 


sion of the UV background of Haardt V Madau (2012) 


(HM2012 hereafter) in the RAMSES code. In the optically 
thin approximation at any given time the photoionisation 
and photoheating rates are independent of location (and 
density). Every part of the simulation sees the same UV 
background intensity. The time evolution of the space av¬ 
eraged UV background intensity is calculated by solving a 
global radiative transfer equation with an empirical opacity 
model and a source function based on the observed UV lu¬ 
minosity function of quasar and star forming galaxies. The 
HM2012 background model calculated in this way should 
give a reasonable approximation to the spatially averaged 
photoionisation rates, especially once the Universe is fully 


ionised (see Puchwein et al. 2015 for a recent discussion). 
By construction, it will however not be able to describe the 
growth of individual HII regions before overlap and the per¬ 
sistence of spatial amplitude fluctuation in the UV flux and 
therefore photoionisation rate for some time after overlap. 

We will later compare such optically thin simulations 
to our full cosmological radiative transfer simulation with 
ATON where the ionising radiation is propagated from dis¬ 
crete sources into the surrounding IGM. Note that we do 
not try to model the ionising radiation from star formation 


in our simulation self-consistently (other than, e.g. Bauer 
et al.|2015| l, but that we assume an ionising volume emissiv- 
ity scaled to that assumed in the HM2012 UV background 
model and identify the dark matter haloes in our simulations 
as the ionising sources. The dark matter haloes are identi¬ 
fied with the HOP halo finder ( Eisenstein V Hut|1998| with 
a minimal halo mass consisting of 10 particles, correspond¬ 
ing to the minimal masses reported in Tables [T| and [2] We 
assume that dark matter haloes above our mass thresholds 
act as ionising sources with emissivities proportional to halo 
mass M similar to the assumptions in |Iliev et al. (2006]), 


N 7 (z) = a(z)M. 


(i) 


Note here that simple models where the UV luminosity 
is proportional to halo mass appear to fit (UV) luminos¬ 
ity functions of high-redshift galaxies remarkably well (see 
Trent S~et al.|[20 10). We will investigate this further in Ap¬ 
pendix [b] 


2.3.2 Choosing the ionising emissivity for the ATON 
simulations 

In order to obtain a realistic post-overlap Lya forest and for 
an easier comparison with our optically thin simulations we 
have scaled our emissivities in the ATON simulations such 
that the volume average are close to that of the HM2012 
model, shown as the black solid curve in Fig. [l] The required 
re-scaling will be discussed in detail below. 

Note again that we thereby have not tried to match 
the space density of sources in the HM2012 model in our 
ATON simulations, but just the redshift evolution of the 
integrated emissivity. For a given population of DM haloes 
identified as ionising sources this fixes then the evolution of 
the efficiency parameter a(z) of equation [l] The HM2012 
model uses a simple escape fraction to model the fraction of 
ionising photons escaping galaxies, 

(f esc ) = 1.8 X 10“ 4 (l + «) 3 ' 4 . (2) 

As ATON is modelling the escape of photons from the 
source positions explicitly we have to account for this and 
boost the escaping ionising emissivity compared to “op¬ 
tically thin” simulations with a spatial homogeneous UV 
background model like HM2012 by a resolution dependent 
factor that accounts for the recombinations in the host 
haloes of our ionising sources. As we will see later “recombi¬ 
nation” boost factors in the range ~ 1.05 — 2 gave reasonable 
reionisation histories. Note that these boost factors should 
be equal to the inverse of the escape fraction in the HM2012 
background model if our simulations were able to correctly 
model the escape of ionising photons from real galaxies. That 
these factors are significantly smaller just means that they 
are clearly and not surprisingly not able to do this. As we 
will see later the redshift evolution of the ionising emissiv¬ 
ity of HM2012 does not get the rather flat evolution of the 
photoionisation rate at 6 > z > 2 as probed by Lya forest 
data quite right. We were able to reproduce the Lya for¬ 
est data and the inferred photoionisation rates by scaling 
the HM2012 ionising emissivities with a redshift dependent 
scaling factor as follows, 
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b(z) = 


a — const 

»[*]" 

[*r [ 


if z ^ z\ 

if Z 2 < z ^ Zi 

if z < Z 2 


(3) 


with suitable values of a, zi, Z 2 , ol i and a 2 as shown in Ta¬ 
bles [l] and [2] Unfortunately, there is no easy way round the 
resolution dependence of the number of recombinations oc¬ 
curring in the source haloes. It required some experimenting 
with the evolution of the ionising emissivity to obtain rea¬ 
sonable reionisation histories for which hydrogen reionisa¬ 
tion is completed at about z ~ 6 and the evolution of r(z) 
is close to the observations at redshifts below six. As the 
monochromatic GPU accelerated ATON radiative transfer 
code is gratifyingly fast (a speedup of 80-100 times com¬ 
pared to CPU-based implementation) we were able to do 
this by trial and error. 


2.4 The simulation set 


We have run a suite of simulations in order to explore the 
convergence in terms of resolution and box size. We also var¬ 
ied the evolution of the ionising emissivity to study the effect 
of the timing of the overlap of HII regions on the Lyo opac¬ 
ity PDF and the spatial fluctuation of the photoionisation 
rate in the post-overlap phase. 

Our standard simulations are based on 512 3 RAMSES 
simulations and we have varied the box size of the simulation 
from 10-40 h _ 1 Mpc. We have also run a 256 3 simulation to 
study the effect of resolution and box-size. 

When varying the reionisation history we have chosen 
a default model which reproduces the measured photoion¬ 
isation rates at z ~ 6 reasonably well, the ‘512-20-good_l’ 
model, where 512 stands for our default number of RAM¬ 
SES resolution elements in one dimension and 20 for the 20 
h -1 Mpc default comoving size of the simulation box. We 
have further run simulations where hydrogen is reionised 
later (model date’) and earlier (model ‘early’) and a model 
where overlap occurs at the same redshift as in the ‘512-20- 
good_l’ model but with a somewhat larger ionising emissiv¬ 
ity (and thus photoionisation rate) in the phase immediately 
past overlap (model ‘good_h’). In tables [l] and [ 2 ] we give the 
basic parameters of our RAMSES simulations as well as the 
parameters of equation [3] Fig. [l] shows the redshift evolu¬ 
tion of the ionising emissivities used in our simulation suite. 
Fig. HI shows the corresponding integrated number of ionis¬ 
ing photons per hydrogen atom. 

Unfortunately, the temperature calculations as imple¬ 
mented in ATON simulations do not properly account for 
the radiative transfer effects on the temperature as the ra¬ 
diative transfer is modelled monochromatically and not fully 
coupled to the hydro simulation. As discussed in detail in 
Puchwein et al. ( 2015| ), (optically thin) simulations with the 
HM2012 UV background model and an equilibrium chem¬ 
istry solver agree actually reasonably well with measured 
temperatures of the IGM. We have therefore assumed the 
temperatures obtained in the RAMSES simulations for the 
ATON simulations. In Fig. [3] we compare the evolution of 
the temperature at mean density as a function of redshift in 


our RAMSES/ATON simulations without radiative transfer 
to measured temperatures from Becker et al. (2011) assum¬ 
ing a power law index 7 = 1.5 for the temperature-density 
relation of the IGM at the relevant densities. The agreement 
is very reasonable. 


3 RESULTS 

In this section we describe the main results from our ATON 
simulations. We first focus on the global reionisation his¬ 
tory of the different models and compare photoionisation 
rates and neutral hydrogen fractions to those inferred from 
Ly a forest data. We then produce mock Ly<a forest spectra 
from our simulations to compare directly to observed prop¬ 
erties of the Ly a forest. We finally make some predictions of 
the UV background fluctuations expected in the later stages 
of the reionisation from our simulations. 


3.1 Reionisation history 

3.1.1 Evolution of the neutral hydrogen fraction 

Figure [4] shows the evolution of the spatial distribution of 
the neutral fraction for the ‘512-20-goodT model. Individ¬ 
ual HII regions start to appear as early as z ~ 16 and by 
z ~ 10 there is considerable overlap of HII regions. Con¬ 
tiguous regions of still neutral gas in the post-overlap phase 
persist to about z ~ 7. Later on fully neutral gas is ex¬ 
pected to be mostly confined to virialized haloes, however 
the ionising sources in our simulations appear to keep the 
gas in their host haloes fully ionised. Note that this may be 
an artefact as the simulation still don’t have the resolution 
to properly resolve the ISM in these haloes. We also do not 
include dust and assume that the haloes continuously emit 
ionising radiation. Note also that we do not have the reso¬ 
lution to properly resolve the majority of mini-haloes with 
virial temperatures below the atomic cooling threshold. The 
smallest haloes in our 512-20 simulation have circular veloci¬ 
ties of I’circ ~ 10.7 km/s somewhat above the atomic cooling 
threshold, but even at this resolution the simulations are 
already significantly incomplete for these haloes (see Ap¬ 
pendix [A] for halo mass functions for our simulation suite). 
Finally note again that the effect of the reionisation heat¬ 
ing is only taken into account by the equilibrium modelling 
with the HM 2012 UV background model in the RAMSES 
simulations without radiative transfer. 

In Fig. [5] we compare the spatial distribution of ionised 
and neutral regions at z ~ 7.1 of our ATON simulations 
(left column) for a range of resolutions and box sizes. The 
ionisation maps are resolution dependent and in particular 
the remaining neutral gas in the haloes in already ionised re¬ 
gions is increasing with increasing resolution as the gas in the 
individual dark matter haloes becomes better resolved and 
is able to recombine and self-shield due to higher densities. 
We show the corresponding RAMSES simulations without 
radiative transfer in the middle column. To make this a fair 
comparison we have rescaled the photoionisation rate in the 
RAMSES simulation to be the same as the mean photoion¬ 
isation rate in the ionised regions in the slice of the ATON 
simulation shown in the left column. 

We have furthermore implemented the self-shielding 











Calibrating cosmological radiative transfer simulations with Lya forest data 5 




(a) (b) 

Figure 1. (a): Evolution of the comoving emissivity of ionising photons (in s -1 Mpc -3 ). The black solid curve represents the evolution 
of the emissivity from HM2012 integrated over all frequencies of ionising photons. The other curves shows the ionising emissivities for 
the different ionisation histories studied with our ATON simulations. The emissivities have been re-scaled relative to the HM2012 model 
to reproduce observed photoionisation rates and to account for the resolution dependent recombinations in the host haloes of ionising 
sources, (b): same as (a) but for the simulations with different box sizes and resolution as described in the text. 



redshift z 



redshift z 


(a) 


(b) 


Figure 2. (a) The cumulative number of ionising photons emitted per hydrogen atoms as a function of redshift for the ionisation histories 
in the ATON simulations, (b): same as (a) but for the simulations with different box sizes and resolution as described in the text. The 
three solid horizontal lines show one, two and three ionising photons emitted per hydrogen atoms at z ~ 6. 


correction suggested by Rahmati et al. (2013) based on their 
radiative transfer simulations in the RAMSES simulations 
without radiative transfer (shown in the right column). In 
the Rahmati et al. (2013) prescription the self-shielding of 
the gas is taken into account by a density-dependent pho¬ 
toionisation rate, obtained by an empirical fit to their radia¬ 
tive transfer simulations (jRah mati et al.| 2013), 


hi = 0.98 x 


+ 0.02 x 



- 2.28 


(4) 


where r ss is the value of the photoionisation rate in 
a cell of the computational box after accounting for self¬ 
shielding while r is the value of the background photoion¬ 
isation rate before accounting for self-shielding. A ss is the 
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model 

512-20-good_l 

512-20-good_h 

512-20-early 

512-20-late 

resolution 

512 3 

512 3 

512 3 

512 3 

box size (Mpc/h) 

20 

20 

20 

20 

minimum halo mass (Mq) 

1.2 x 10 8 

1.2 x 10 s 

1.2 x 10 s 

1.2 x 10 s 

a 

1.3 

1.3 

1.4 

1.2 

z\ 

5.8 

7 

5.8 

5.5 

Z2 

5 

5 

5 

5 

ai 

0.6 

0.5 

0.6 

0.9 

«2 

1.0 

1.1 

1.0 

0.8 


Table 1. Basic properties of the simulations for the 512-20 ATON simulations with the different reionisation histories. 


model 

256-10 

512-40 

512-10 

resolution 

256 3 

512 3 

512 3 

box size (Mpc/h) 

10 

40 

10 

minimum halo mass (Mq) 

1.9 x 10 8 

8.9 x 10 8 

1.6 x 10 7 

a 

1.36 

1.025 

1.58 

zi 

6 

5.8 

6 

Z2 

5 

5 

5 


0.75 

0.5 

1.0 

OL2 

0.85 

1.0 

1.0 


Table 2. Basic properties of the simulations for the simulations with different box sizes and resolution. 


overdensity above which the gas begins to self- shield (see 
Bolton &; Haehnelt]|2013 ), 


(WM 1/3 1 

( fe ) 

t - 2/3 /l + *\- s 

V0.61/ ' 

U-08J 

! V 8 J 


A ss = 36 x T% 3 T? /16 


_ (5) 

where, A = p/p is the overdensity, Ti 2 = r/10“ 12 s _1 is 
the background photoionisation rate, T 4 = T/10 4 K, with T 
the temperature of the gas, fi is the mean molecular weight 
and f e = n e /riH is the free electron fraction with respect to 
hydrogen. 

We have chosen z — 7.1 for this comparison as there is 
particular interest in this redshift due to the highest redshift 
QSO ULAS J1120+0641 at z — 7.085 which may show signs 


of a red damping wing in its absorption spectra (Mortlock 
|et al.| 20 n |Bolton et al.| 20 lH but see Bosman et al. 2015 ( 
prep)). Note that while there are still substantial contiguous 
neutral regions extending well beyond the extent of individ¬ 


ual haloes in the ATON simulations the volume filling factor 
of the neutral gas in the optically thin simulations is con¬ 
siderably smaller and is confined to the haloes even when 
self-shielding with the Rahmati et al. (2013) prescription is 
taken into account. The self-shielding of dense gas in haloes 
located in already ionised regions of the ATON simulations 
is considerably weaker than predicted by the | Rahmati et al.| 
(2013) model which does not account for the local ionising 
sources located in these haloes which are emitting with con¬ 
stant emissivity in our models. 


Figure [ 6 ] shows the evolution of the spatial distribution 
of the neutral fraction at three redshifts for the four ‘512- 
20 ’ models with the corresponding emissivities reported in 
Table ^ As expected, increasing/decreasing the emissivities 
accelerates/delays reionisation. The different models bracket 
a reasonable range of reionisation redshifts. The differences 
in the timing of the reionisation process between the differ¬ 
ent models are further quantified in the upper left panel of 
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Figure 3. Evolution of the temperature at mean density To in 
our different RAMSES/ATON simulations. The data points show 
observational constraints from [Becker et al.| (20111) assuming a 
value of 7 = 1.5 for the slope of the temperature-density relation. 


Figure [7] that shows the evolution of the mean neutral frac¬ 
tion, while the bottom left panel of the same figure shows 
the evolution of the volume filling factor for our suite of 
simulations. As already mentioned our ‘good’ models were 
chosen such that hydrogen reionisation is mostly completed 
by z = 6 — 6.5. At this redshift the still neutral regions cease 
to hang together and become separate islands and the mean 
neutral fraction drops rapidly by several orders of magni¬ 
tudes to a level of /hi ~ 10 -4 , while the volume filling factor 
of ionised regions rises to Qhii = 1. The same happens at 
somewhat lower/higher redshift in our late/early model. 

Figure [8] shows the integrated Thompson optical depth 
for our simulation suite, 


t(z) = ca t 



(6) 


where at is the Thompson cross section of the electron, 
and n e (z) is the electron density. At low redshift, where the 
Universe is fully ionised the curves all converge. Our reion¬ 
isation histories are all somewhat below the mean value of 
the 2013 Planck+WMAP data results and are in excellent 
agreement with the Planck 2015 TT+lowP+lensing results. 
There is little dependence of the Thomson optical depth 
on resolution as we have adjusted the emissivity to obtain 
similar reionisation redshifts. The slightly earlier reionisa¬ 
tion in the ‘256-10’ model results in a somewhat larger op¬ 
tical depth as expected. Note that we have assumed here 
that He becomes singly ionised in regions where hydrogen 
is ionised. Note further that the reionisation history at high 
redshift (say z > 10) is poorly constrained by our compari¬ 
son to Lya forest data. There should thus exist models with 
a different evolution of the ionising emissivity at very high 
redshift (where we have not modified here the redshift evo¬ 
lution assumed in HM2012) that are equally consistent with 
the Lya forest data and have a different Thomson optical 
depth. 


3.1.2 Photoionisation rate and mean free path of ionising 
photons 

The coloured curves in Fig.[9]show the evolution of the space 
averaged hydrogen photoionisation rate for our ATON sim¬ 
ulations. Note that we have excluded not yet ionised re¬ 
gions when calculating the average photoionisation rates. 
The photoionisation rate is first decreasing as the HII re¬ 
gions expand and the distances from the ionising sources 
increases and it then increases once overlap is achieved and 
multiple sources contribute to the local photoionisation rate. 
By moderately scaling our emissivities with regard to the 
HM2012 UV background model we succeeded in reproduc¬ 
ing the observed photoionisation rates rather well. The only 
exception is the ‘512-20-late’ model that has photoionisation 
rates at z ^ 5 somewhat lower than observed. Reionisation 
appears indeed to happen too late in this model to be con¬ 
sistent with the observed photoionisation rates as we had 
intended for this model. Tuning the ionising emissivities to 
reproduce the observed apparent lack of evolution of the 
photoionisation rate at 6 > z > 2 in the other models re¬ 
quired a significantly shallower rise of the ionising emissivity 
compared to the HM2012 model at 5 > z > 2 (see Fig. 1) 
that assumes a rather rapid rise of the ionising emissivity due 
to QSOs in this redshift range (see Tables [l] and [5] for the co¬ 
efficients of the power law of equation [3]) . Note that Munoz 
et al. (2014) came to similar conclusions based on their an¬ 


alytic modelling of photoionisation rates. Our ‘good_l’ and 
‘good_h’ models nicely bracket the inferred photoionisation 
rate at 6 > z > 5 where the photoionisation rate is most 
difficult to measure and also evolves rather rapidly. 

To better understand the evolution of the photoionisa¬ 
tion rates it is illustrative to also have a look at the evo¬ 
lution of the mean free path of ionising photons relative to 
the mean distance between its sources which we show in 
Fig. |10| In all models the mean free path shown by the solid 
coloured curves is initially shorter than the mean distance 
between ionising sources and rises slowly with decreasing 
redshift until mean free path and mean distance between 
sources become similar. Soon afterwards the ionised regions 
fully percolate and the mean free path rises very rapidly. The 
evolution of the photoionisation rate then scales as the prod¬ 
uct of ionising emissivity and mean free path, L oc eionA m f p . 
The photoionisation rate therefore shows a similar rapid rise 
when the ionised regions fully percolate. This behaviour of 
r and A m f p has already been discussed for the early cos¬ 
mological radiative transfer simulations by Gnedin (2000). 
There is good agreement with observations of the mean path 


as compiled by Worseck et al. (2014), but note that as dis¬ 


cussed by Becker V Bolton (2013) at z < 4 the approxima¬ 


tion made here that the mean-free path is small compared to 
the Hubble radius becomes rather poor. At z < 4 redshifting 
of ionising photons between emission and absorption has to 
be taken into account for a proper comparison of the mean 
free path predicted by the simulation and that inferred from 
Lya forest data. 

For reference we also show the mean free path corre¬ 
sponding to the assumed opacity in the HM2012 UV back¬ 
ground model. There is a clear trend of increasing mean free 
path (as well as photoionisation rate) with increasing box 
size and decreasing resolution. Full details about how the 
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Figure 4. Evolution of the spatial distribution of the neutral fraction for the ‘512-20-good_l’ model for a slice of 39.06 comoving kpc 
thickness in the mid plane of the simulation. The white pentagons indicate the location of the dark matter haloes assumed as ionising 
sources with the size scaled to the mass of the halo. 
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Figure 5. Spatial distribution of the neutral fraction at z = 7.1 : First, second, third and fourth rows are for the ‘256-10’ simulation, 
the ‘512-40’ simulation, the ‘512-20-good_l’ simulation and ‘512-10’ simulation, respectively. The first column shows the full radiative 
transfer simulations with ATON. The second column shows the optically thin RAMSES simulations with photoionisation rate rescaled 
to match that of the corresponding ATON simulation in the first column. The third column shows the RAMSES outputs post-processed 
with the model/fit of |Rahmati et al. |2013] > to account for self-shielding to ionising photons. The photoionisation rate in the slices shown 
is indicated on the plot. 
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Figure 6. The spatial distribution of the neutral fraction at z = (7.5, 6.3, 5.7) for the different ionisation histories in our ATON 
simulations. First, second and third rows corresponds to the ‘512-20-early’, ‘512-20-good_r and the ‘512-20-late’ models. The white 
pentagons indicates the location of the dark matter haloes assumed as ionising sources with the size scaled to the mass of the halo. 


mean free path is measured from the simulations are given 
in Appendix [C] 

3.2 Spatial UV background fluctuations in the full 
radiative transfer simulations with ATON 

Here we will have a closer look at the spatial fluctuations of 
the photoionisation rate in our full radiative transfer simu¬ 
lations. In Fig. mi we show the PDF of the phot ionisation 
rate F ±2 in our different ATON simulations for a range of 


redshifts as well as the value of F \2 of the HM2012 model 
used in the RAMSES simulations without radiative transfer 
at those redshifts. 

The first thing to note is that, at low redshift, the PDF 
of the photoionisation rate is strongly peaked and that the 
spatial fluctuations are rather small. With increasing red- 
shift the distributions broadens considerably and a signifi¬ 
cant tail towards low values of F develops. 

As expected the damping of the spatial fluctuations cor¬ 
relates with the timing of reionisation and (moderate) spa- 
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(a) (b) 




redshift z redshift z 


(c) (d) 

Figure 7. (a) Evolution of the mean volume-weighted neutral fraction for the ATON simulations with different reionisation histories, 
(b) same as (a) for the simulations with different box sizes and resolutions, (c) evolution of the volume filling factor of ionised regions 
(with an ionisation fraction ^ 0.5) for the ATON simulations with different reionisation histories, (d) same as (c) for the simulations 
with different box sizes and resolutions. The data points are from|Fan et al. (2006| and are based on an analytical model for the opacity 
PDF. 


tial fluctuations persist for some time after full percolation 
has occured. 

The trend of increasing mean photoionisation rate with 
increasing box size and decreasing resolution is again clearly 
recognisable. Other possible trends in the PDF of T with the 
simulations of different box size and resolution appear to be 
weak. 


3.3 Reproducing Lyman-a forest data 

3.3.1 Extracting mock Lyman-a forest spectra 

We have produced mock absorption spectra from our simula¬ 
tions for comparison with Lya forest data as e.g. described 
in |Theuns et al.|[T998| For this we used the gas density, 


gas temperature, gas velocity and ionisation state of the gas 
in the simulation boxes and compute Lyman-a absorption 
spectra along random lines of sight. We did this for both the 
RAMSES simulations without the radiative transfer and the 
ATON simulations. Note again that the temperatures of the 
(optically thin) RAMSES simulation were used throughout 
as discussed in section [274] We have extracted 5000 spectra 
from each simulation output. For comparison with the ob¬ 
served Lya opacity PDF we concatenated simulated spectra 
to a length covering a comoving distance of 50 Mpc/h, the 
length chosen by Becker et al. (2015) for their measurements. 
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(a) (b) 

Figure 8. (a) The integrated Thompson optical depth for the ATON simulations with the different reionisation histories, (b): Same as 
(a) for the simulations with different box sizes and resolution. The solid black line with the light blue area represents the 68% confidence 
limits from the 2013 Planck+WMAP data (see |Planck C ollabora tion et al.] 2014). The dashed black line with the shaded dark blue area 
represents the 68% confidence limits from the Planck 2015 TT+lowP+lensing data release (see |Plan ck Collaborat ion et al.|2015| >. 




(a) (b) 

Figure 9. (a) Evolution of the hydrogen photoionisation rate for the ATON simulations with different reionisation histories (computed 
in ionised regions with an ionisation fraction ^ 0.5). (b) same as (a) for the simulations with different box sizes and resolutions. The solid 
black line shows the evolution of the photoionisation rate in the uniform HM2012 UV background model. The observational constraints 
are from |Bolton &; Haehnelt| ( |2007| >, |Calverley et ah] ( |2011| ), |Wyithe &; Boltonj ( |2011]) and |Becker &; Boltonj |2013[ ). Note that the data 
points have been rescaled (by factors ranging from 0.85 to 1.03) to that inferred for the cosmological parameters and temperature-density 
relation in our RAMSES simulations with the scaling relation in [Bolton et al. (2005J and Bolt on Sz Haehnelt| |2007|. 


3.3.2 Effective optical depth 

We expect our simulations to reproduce the properties of 
the 2 < z < 4 Ly a forest data reasonable well as our op¬ 
tically thin 512 3 RAMSES simulations are very similar to 
P-GADGET3 simulations that some of us have compared 
extensively to Lya forest data in the past. 

In Fig. [T2] we compare the evolution of the effective op¬ 


tical depth T e ff(z) calculated from 5000 spectra for each of 
our simulations with a range of measurements of the effec¬ 
tive optical depth as described in the figure caption. Both 
‘good’ models reproduce the rapid rise of the observed opti¬ 
cal depth at z ~ 5.5 very well, while in the ’early’ and ’late’ 
models reionisation is completed indeed somewhat early and 
late, respectively. The effective optical depth of our ‘good_h’ 
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redshift z 



(a) 


(b) 


Figure 10. (a) Evolution of the space averaged mean free path of Lyman limit (912 A) photons (in proper Mpc) for the ATON 

simulations with different reionisation histories (including both ionised and neutral regions), (b) same as (a) for the simulations with 
different box sizes and resolutions. The dashed curves show the evolution of the mean distance between ionising sources in the different 
ATON simulations. The black squares with errorbars show observational constraints from O’ Meara et al.| ( |2013| ) at z = 2.44, |Fumagalli| 
|et a.l.| < [2013^ at z = 3.00, |Prochaska et al.| ( |20d9] > at z = 3.73, and |Worseck et ah] ( |2014| > at z = 4.56 as compiled by |Worseck et al.| ( |2014| ). 
The black solid line shows the evolution in the HM2012 model. 


model goes thereby right through the middle of the observed 
measurements of r e ff which however show a large scatter. 
The effective optical depth of the ‘good_l 5 model on the 
other hand goes through the upper envelope of observed 
measurements at 5 < z < 6 as expected from our compar¬ 
ison of simulated photoionisation rates and those inferred 
from Lya forest data in section |3.1.2| Note that our full 
radiative transfer simulations with the rescaled emissivities 
compared to the HM2012 reproduce the observed r e ff signif¬ 
icantly better than the optically thin RAMSES simulations 
with the HM2012 model. The latter somewhat overproduces 
the observed r e ff at 4 < z < 6 and shows a somewhat slower 
evolution with redshift than observed. This confirms similar 


findings by Puchwein et al. (2015) based on optically thin 


P-GADGET3 simulations. The ionising emissivities in the 
HM2012 UV background model at 5 < z < 6 appear to be 
somewhat lower than necessary to reproduce the Lya for¬ 
est data. Judging from Fig. [l] the rather rapid increase of 
the ionising emissivity at 3 < z < 6 due to QSOs appears 
at fault here. An increased contribution from faint AGN at 
z > 4 (e.g. Giallo ngo et al.||2015 ) could provide a redshift 


evolution in better agreement with the opacity data. 


3.3.3 The incidence of Lyman limit systems 

Another important test of the (residual) neutral hydrogen 
distribution regulating the Lya opacity are Lyman Limit 
Systems in QSO absorption spectra. For this purpose, we 
calculated the evolution of the incidence rate of absorbers 
dN/dX where X is the usual absorption distance, 


X = 


(7) 


and H(z) is the Hubble constant at redshift z (Bahcall X 
Peebles|| 1969 ). 


In practice we search in our mock spectra for all op¬ 
tically thick pixels (r > 3). We compute then the column 
density of each connected absorber defined in that way and 
select those with N hi ^ 10 17 cm -2 as Lyman Limit Sys¬ 
tems (LLSs). We thereby exclude regions not yet ionised 
where the occurrence of LLSs is poorly defined. 


In Fig. 13 we compare the evolution of dN/dX for all 
our simulations to the observed incidence rate from a cumu¬ 
lation of Lya forest data as described in the figure caption. 
The solid curves show the results for the ATON simula¬ 
tions while the dashed and dotted curves are respectively 
for the optically thin simulations without and with the Rah- 


mati et al. (2013) self-shielding correction (see also Kohler 
X Gnedin) 2007 and McQuinn et al.||2011 ) 


Overall the simulations reproduce the increase of the in¬ 
cidence rate of LLSs with increasing redshift very well, but 
underpredict the observed incidence rate by about a factor 
1.5-2 for the ATON simulations as well as for the RAM¬ 
SES simulations without radiative transfer. This is similar 
to what has been found by other studies in the literature. 
In optically thin simulations without (supernovae) feedback 
the neutral hydrogen distribution in galactic haloes is gen¬ 
erally found to be too spatially concentrated to reproduce 
the observed LLS incidence rate ( Schaye et al.||2010 Altay 
et al.|2013 Bird et al.|2014 ). 


Fumagalli et al. (2011) used high-resolution zoom sim¬ 


ulations to demonstrate that local stellar radiation further 
reduces the abundance of high column density absorbers. 
Similar conclusions have been reached by |Rahmati et al.| 
(2013) who found that radiative transfer simulations with 


local stellar sources underpredict the rates of incidence of 
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Figure 11. Probability distribution function (PDF) of the photoionisation rate (dP/dlog 10 ri 2 ) at five different redshifts in our ATON 
simulations with different reionisation histories (left column) and for the simulations with different box sizes and resolutions (right 
column). The solid black vertical lines show the values of Ti 2 in our RAMSES simulations without radiative transfer based on the 
HM2012 UV background model (their Table 3) at the same redshifts. The dashed black curve in the bottom left panel shows the PDF 
of our ‘bright source’ model discussed in section [3. 3. 4| and Appendix [P] 


absorbers with column densities 10 19 < Nhi < 10 21 cm 2 
by factors of a few. 

Simulations with the same resolution but different box 
size (i.e. the 512-20-good and the 256-10 models) show very 
similar evolution. There is, however, a clear increase of the 
incidence rate of LLSs with increasing resolution at z > 2.5 
and our highest resolution 512-10 model starts to match the 
data reasonably well. Part of the discrepancy between ob¬ 


served and simulated LLS incidence rate is therefore clearly 
due to insufficient resolution and not the lack of stellar feed¬ 
back in our simulations. Note that our RAMSES simulation 
with the Rahmati self-shielding prescription tend to give 
slightly larger LLS incidence rates than our full radiative 
transfer simulations with ATON where ionising sources are 
located within the LLS. 
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Figure 12. ( (a) The coloured solid curves show the evolution of the of effective Lya optical depth r e ff(z) for the ATON simulations 
with the different ionisation histories, (b) same as (a) for the simulations with different box sizes and resolutions. The dashed curves 
corresponds to the evolution of r e ff (z) for the RAMSES simulations without radiative transfer. The black data points show the observed 
r e ff for the QSOs presented in |Fan et al.|2006] the green data points show observational constraints from |Becker et ah] ( | 201 3) and the 
blue data points show the observed r e ff for the QSOs presented in |Becker et al.| ( |2015| ). 




(a) (b) 

Figure 13. (a) The coloured solid curves show the evolution of the incidence rate dN/dX of LLSs for the ATON simulations with 
the different ionisation histories, (b) same as (a) for the simulations with different box sizes and resolutions. The dashed curves are for 
the RAMSES simulations and the dotted curves are for the RAMSES simulations post-processed with the self-shielding correction of 
|Rahmati et al.| ( |2013| >. For the RAMSES and RAMSES with self-shielding correction simulations, the neutral fraction of the gas has been 
rescaled to be consistent with the mean T value of the corresponding ATON simulations. The black data points with errorbars show the 
measured incidence rate of LLSs (from [Son gaila Sz Cowie|201 0 and Fumagalli et al. 2013| >. 


3.3.4 The cumulative Lya opacity PDF 


Recently, Becker et al. (2015) have presented accurate mea¬ 


surements of the PDF of the Lya effective optical depth 
obtained from chunks of spectra of 50 comoving Mpc/h at 
4 < z < 6 based on a large sample of high-quality high- 
redshift QSO absorption spectra. Based on a comparison 


with P-GADGET3 simulations without radiative transfer 


similar to our RAMSES simulations Becker et al. 


(2015) 


argued that at z > 5 the Lya opacity PDF deviates increas¬ 
ingly from that expected for a spatially homogeneous UV 
background and attributed the increasing deviation at the 
high opacity end to the spatial UV fluctuations expected 
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T eff T eff 


Figure 14. The coloured solid curves show the cumulative PDF of the effective optical depth r e ff for our ATON simulations at ten 
different redshifts as indicated on the plot. The dashed curves show the corresponding distributions for the RAMSES simulations without 
radiative transfer. The a a and <ar values in the labels of the different panels are the factors by which we have rescaled the optical depth of 
the ATON and RAMSES simulations to match the measurements of |Becker et al.| p015) shown by the solid black curves at P(r e ff) = 0.5. 
The dotted curves show the corresponding distributions of the ATON simulations before the rescaling. The two bottom panels show 30 
individual realisations (solid yellow curves) of the cumulative flux PDF for the two highest redshift bins for samples of simulated spectra 
from the 512-20-good_l ATON simulation with total path length equal to that of the observed sample. 
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Figure 15. Same as figure [l4| for the four highest redshift bins for simulation of different box size and resolution. 



Figure 16. Upper panel : the coloured solid curves show the cumulative PDF of the effective optical depth r e ff for our ATON simulations 
for the two highest redshifts as indicated on the plot. The effective optical depth is now computed for chunks of a length corresponding 
to the box size in the different simulations instead of chunks of length 50 Mpc/h as before. The two bottom panels show 30 individual 
realisations (solid yellow curves) of the cumulative flux PDF for the two highest redshift bins for samples of simulated spectra from the 
’512-10’ ATON simulation with total path length equal to that of the observed sample. The effective optical depth is again computed 
for chunks equal to the box size of the simulation, 10 Mpc/h in that particular case. 




















































18 J. Chardin et al. 


Becker 2015 
RAMSES 
a = 0.47 

ATON '512-20-good_h 
a = 0.80 
Bright source 

A mfp =35.0 Mpc/h 
N bright = 2 
/3 ga i = l-15 
bright=1-25 



Figure IT. The red solid curve shows the cumulative effective op¬ 
tical depth PDF for a “bright source model” with two bright ion¬ 
ising sources in a periodically replicated volume of (100Mpc/h) 3 
contributing significantly to the ionising UV background at red- 
shift z ~ 5.8. The solid yellow curves show 30 individual realisa¬ 
tions of the cumulative flux PDF for samples of simulated spec¬ 
tra with photoionisation rates of the toy model with total path 
length equal to that of the observed sample. The assumed con¬ 
tribution from rare bright sources is added to that in our 512-20- 
good.l ATON simulation such that the mean photoionisation rate 
is equal to that in the corresponding good.h model. The photoion¬ 
isation rates are then scaled by a global factor 0.87 to match the 
observed PDF as measured by |Becker et al.| (2015). An average 
mean free path of ~ 35 comoving Mpc/h as in our 512-20-good_h 
ATON model is assumed for the attenuation of the bright sources. 
The black dotted and dashed curves show the PDF of the 512- 
20-good_h ATON simulation and the PDF of our corresponding 
optically thin RAMSES simulation with optical depths rescaled 
to match the observed PDF at P(< r e ff) = 0.5. For more details 
see text and appe ndix [P] The black solid curve shows the data 
from |Becker et a l. ( |2015| >. 


to persist for some time during the post-overlap phase of 
hydrogen reionisation. As we have discussed in section 3.2 
our full radiative transfer simulations with ATON show such 
spatial fluctuations, albeit at a rather moderate level. 

In Fig. [T4| and [15] we compare the cumulative Lya opac¬ 
ity PDF at 4 < z < 6 for our ATON and RAMSES simu¬ 
lation without radiative transfer with the measurements of 


Becker et al. (2015). As discussed in section 3.3.2 our ATON 


simulations for the ‘good_h’ model, shown here as the red 
dotted curves, reproduce the mean effective optical depth 
best. At high redshift the observed effective optical depth at 
a given redshift shows, however, clearly a much larger scat¬ 
ter than in our simulations. The discrepancy is similar to 
that found by Becker et al. (2015), for their optically thin 
simulations. 

To investigate this further we adjust the optical depth 
of our simulations by a constant factor such that the mean 
effective optical depth of the simulated spectra is equal to 
r e ff of the observed distribution at P(r e ff) = 0.5. The cor¬ 
rection factor a by which the optical depth is multiplied is 
indicated for each simulation on the plots. Let us now again 


have a look at Fig. |14| where we compare the different reion¬ 
isation histories. At z < 4.8 the Ly<a opacity PDF from the 
ATON and optically thin RAMSES simulations are almost 
identical after rescaling. They also agree very well with the 
observed PDF, very similar to what was found by |Becker| 


et al. (2015) based on their P-GADGET3 simulations. At 
z ~ 5 and z ~ 5.2 the PDF of the simulated spectra be¬ 
gins to exhibit a somewhat steeper slope than the observed 
PDF. At z > 5.4 the discrepancy increases and especially in 
the two highest redshift bins observed and simulated PDF 
of the effective optical depth are strikingly different. 

The observed sample is still rather small. In order to 
get a feel for the expected sample variations the two bottom 
panels of Fig. [14] show, for the two highest redshift bins, 
the observed Lya opacity PDF with 30 individual realisa¬ 
tions for samples of simulated spectra of the ‘520-20-goodJi’ 
model with total path length equal to that of the observed 
sample. The scatter is substantial, but the discrepancy be¬ 
tween observed and simulated PDF is clearly much larger 
than can be explained with sample to sample variations. 

In Fig. [l5j we show the Ly<a opacity PDF for simu¬ 
lations where we have varied box size and resolution. The 
differences are moderate. As expected increasing the box 
size somewhat increases the width of the distribution, but 
not nearly enough to attribute the discrepancy with the ob¬ 
served PDF to an insufficient box size of our simulations. 

That the PDF of the effective optical depth of our 
ATON and RAMSES simulations are rather similar means 
that the spatial UV fluctuations in our full radiative sim¬ 
ulations are not sufficient to explain the observed rather 
wide distribution of optical depth at 5 < z < 6. |Becker| 


et al. (2015) modelled the expected UV fluctuations in the 


post-overlap phase for ionising emission from (faint) galax¬ 
ies. Combining a simple model for the attenuation of ionising 
photons for plausible assumptions for the mean free path 
with their optically thin simulations they came to similar 
conclusions. 


Becker et al. (2015) therefore proposed that the increas¬ 


ing tail of large optical depth in their data could be due to 
large spatial fluctuations of the mean free path of ionising 
photons expected in the post-overlap phase of reionisation 
(Furlanetto V Oh (2005), Mesinger V Furlanetto 2009). 

To shed further light on this in Fig. [l6j we look at the 
influence of the chunk size used to compute the effective op¬ 
tical depth. We use our simulations of different box sizes for 
this and set the chunk size equal to the size of the simulation 
box. Reducing the length of the spectrum used to compute 
the effective optical depth increases the scatter in r e ff and 
thus the width of the PDF as expected. Interestingly for a 
chunk/simulation box size of 10 Mpc/h the optical depth 
PDF matches the observed PDF (calculated with a chunk 
size corresponding to 50 Mpc/h) to within the expected sam¬ 
ple variation. Our simulations appear to reproduce the ob¬ 
served range of observed optical depth at a given redshift, 
but the coherence length of the optical depth appears to be 
much shorter than observed so that averaging reduces the 
width of the PDF. 

Density fluctuations on these rather large scales are very 
small. However, as we will show below large scale spatial 
UV fluctuations that extend over 50 Mpc/h scales or more 
could explain this large coherence length. As we and [Becker] 


et al. (2015) have shown, such fluctuations are, however, 
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at odds with a UV background at 5 < z < 6 dominated 
by the large number of faint galaxies widely believed to be 
responsible for reionisation. Leaving this point aside here for 
the moment (we will come back to this in the next section) 
we explore here a (toy) model (see also Appendix [d]), where 
the galaxies driving reionisation contribute only a fraction 
of the UV background at 5 < z < 6, while the remainder is 
contributed by much rarer sources (faint AGN or perhaps 
very bright galaxies). 

Figure [lT] shows the resulting effective optical depth 
PDF for a model with two bright sources in a (periodically 
reproduced) cube with side-length 100 Mpc/h. We have used 
the photoionisation rates in our c 512-20-good_l’ model and 
have added the UV emission from the rare bright sources 
such as to reproduce the mean T of the 512-20-good_h ATON 
model which is the closest to the observations. The result¬ 
ing large scale fluctuations of the ionising flux significantly 
widen the optical depth PDF. As Figure [T7| shows this par¬ 
ticular set of parameters reproduces the observed PDF at 
z = 5.8 very well. 

The large scale fluctuations of the ionising flux due 
to the rare sources in this model were calculated with the 
simple attenuation model first introduced by Bolton et al.| 
(2006) and used by Becker et al. (2015) (but for a much 
higher space density than considered here). For details and 
limitations of this toy model we refer the reader once more to 
Appendix[D] When choosing the parameters of the model we 
found that the optical depth PDF depends very sensitively 
on the space density of sources, mean free path of ionis¬ 
ing sources and fraction of the UV background contributed 
by rare sources and there will certainly be other parameter 
choices that will also reproduce the data. We should thus 
emphasise here that our toy model should just be seen as 
proof of concept. We will leave a more detailed exploration 
of the effect of spatial UV fluctuations due to rare bright 
sources on the effective optical depth PDF to a future pa¬ 
per. 


3.4 Implications for the role of bright galaxies 

and (faint) AGN for the origin and evolution 
of the ionising emissivity at high redshift 


than 50Mpc/h. The very extended Gunn-Peterson trough 
of high Lya optical depth in ULAS J0148+0600 at z ^ 5.5 
is particularly striking in this regard. 

Addressing this properly will require much larger simu¬ 
lated regions at similar resolution(s) than we have employed 
here and will be a formidable computational task. As we 
have demonstrated with our toy model in the last section a 
significant contribution to the integrated ionising UV back¬ 
ground from bright sources with space densities of order 
10 -6 (Mpc/h) -3 appears to be required. 

While most of the literature including HM2012 have 
assumed that the contribution of QSOs drops rapidly at 
z > 4, some authors have argued that there may be a signif¬ 
icant contribution by AGN to the ionising UV background 
at 2 > 4 (see Giallongo et al. 2015 for a recent discussion). 
As discussed by Haardt V Salvaterra (2015) a contribution 
of (faint) AGN at z > 4 is consistent with the soft X-ray 
background, but sources with softer spectra would have to 
be responsible for driving hydrogen reionisation in order not 
to exceed measurements of the soft X-ray background. In 
light of the need for spatial fluctuations of the ionising UV 
background on scales > 50 Mpc/h at 5 < z < 6 it appears 
therefore prudent to also consider the possibility of efficient 
escape of ionising photons from stars in rare very bright 
(possibly star bursting) galaxies at z > 5. These galaxies 
would, however, have to have unusually hard ionising spec¬ 
tra and large (of order unity) escape fractions to provide the 
required ionising emissivity. Interestingly the recent discov¬ 
ery of strong CIV emission in high-redhift galaxies by [Stark] 


et al. (2015) may point in this direction. In either case the 


current widely accepted assumption that the ionising UV 
background at z > 5 is dominated by faint sub-L* galaxies 
may need revision. 


4 CONCLUSIONS AND OUTLOOK 

We have calibrated here full cosmological radiative transfer 
simulations performed with ATON/RAMSES with Lya for¬ 
est data and compared them with optically thin RAMSES 
simulations. 

Our main results are as follows. 


As discussed in detail in the previous sections with our 
ATON simulations we have got good agreement with the 
current observational constraints from Lya forest data in 
the redshift range 2 < z < 6. With a moderate adjust¬ 
ment of the evolution of the ionising emissivity compared to 
HM2012 we have been able to reproduce the observed evo¬ 
lution of the mean neutral fraction, mean photoionisation 
rate, averaged mean free path, mean Lya opacity and the 
incidence rate of LLSs. 

Similar to results from optically thin simulations we 
failed, however, with our radiative transfer simulation to re¬ 
produce the broad PDF of the optical depth at 5 < z < 6 
when measured for chunks of 50 Mpc/h as in Becker et al. 


(2015) . We have thereby confirmed the result of Becker et al. 


(2015) that the corresponding optical depth fluctuations on 
these large scales are much larger than can be plausibly ex¬ 
plained as being due to density fluctuations. Judging from 
the large differences in the visual appearance of individual 
QSO spectra, the implied large coherence length of order 
unity optical depth fluctuations appears to be even larger 


• Once the amplitude of ionising emissivities in the 
HM2012 UV background model is rescaled to take into 
account (the resolution dependent) recombination in the 
haloes of ionising sources, we were able to obtain reioni¬ 
sation histories for our full radiative transfer simulations in 
very good agreement with a wide range of available obser¬ 
vational constraints derived from Lya forest data: mean ef¬ 
fective optical depth, mean free path of ionising photons, 
incidence of LLSs, photoionisation rates. 

• To reach good agreement with the rather flat evolu¬ 
tion of inferred photoionisation rates at 2 < z < 6 re¬ 
quired, however, a moderate re-scaling of the evolution of the 
ionising emissivity with redshift compared to the HM2012 
model. The HM2012 model tends to have a too steep red¬ 
shift evolution of the ionising emissivity in the redshift range 
3 < z < 6, while in the same range its normalisation is 
too low to produce the photoionisation rates inferred from 
Lya forest data. 

• In simulations that match the Lya forest data well, per¬ 
colation of individual HII regions occurs at z ~ 7 and mod- 
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erate spatial fluctuations of the photoionisation rate persist 
to about z ~ 6. Reionisation histories where percolation oc¬ 
curs sufficiently late to allow for larger spatial fluctuations of 
the photoionisation rate at z < 6 appear to be inconsistent 
with the Lya forest data. 

• Reproducing the broad PDF of the effective optical 
depth in chunks of 50 comoving Mpc/h at redshift z ~ 5 — 6 
proved rather difficult. Spatial UV fluctuations on much 
larger scales (> 50 comoving Mpc/h) than produced by (sub- 
) L* high-redshift galaxies appear to be required to repro¬ 
duce the observed PDF of the optical depth suggesting a 
contribution of order unity to the ionising UV background 
by sources with space densities of order 10 -6 (Mpc/h) -3 . 
This is very different from the widely accepted assumption 
made in HM2012 that the ionising emissivity at these red- 
shifts is dominated by very faint sources, but is consistent 
with recent suggestions for a significant contribution of AGN 
to the ionising emissivity at z > 4. Such an increased contri¬ 
bution from faint AGN at z > 4 could at the same time pro¬ 
vide a redshift evolution of the ionising emissivity in better 
agreement with the redshift evolution of opacity data. Bet¬ 
ter characterising the coherence length of the spatial fluctu¬ 
ations of the Ly<a optical depth at z > 5 should provide valu¬ 
able further information on the contribution of rare bright 
sources to the ionising emissivity at high redshift. 

Further progress will require full radiative transfer sim¬ 
ulations of much larger regions at a similar resolution as 
that of the 512-20 models discussed here that include rare 
bright galaxies and AGN. Another important prerequisite 
for further progress will be the realistic modelling of the 
angular and temporal distribution of the ionising radiation 
from starbursts and AGN including light travel time effects. 
Such simulations of larger regions should hopefully enable us 
to obtain realistic predictions of the (coupled) spatial fluctu¬ 
ations of mean-free path of ionising photons and photoioni¬ 
sation rates on the relevant large spatial scales. It will also 
be important to further improve the modelling of the gas 
distribution in galactic haloes by including stellar feedback 
in order to better reproduce the incidence of LLSs acting as 
sinks of ionising radiation. 

Other avenues for improvement are multi-frequency 
treatment of the radiative transfer with the aim of an im¬ 
proved modelling of the temperature evolution of the IGM 
during hydrogen and helium reionisation, the self-consistent 
coupling of the effects of radiation pressure and the photo¬ 
heating of the gas and eventually a self-consistent modelling 
of the ionising sources including possible self-regulating ef¬ 
fects due to the reionisation of hydrogen. 

Overall we have shown here that Lya forest data pro¬ 
vides very strong constraints on how reionisation proceeeds 
and we expect Ly<a forest data to become an even more im¬ 
portant calibrator of further improved full radiative transfer 
simulations of the reionisation of hydrogen as we push to 
higher redshift with ongoing and planned QSO surveys. 
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APPENDIX A: THE HALO MASS FUNTION IN 
THE RAMSES SIMULATIONS 


We have identified the dark matter haloes in our RAM¬ 
SES simulations with the sources of ionising photons in our 
ATON simulations. In this appendix, we present the evolu¬ 
tion of the halo mass function of our RAMSES simulations 
with different resolution and box sizes (Fig. Al). For com¬ 
parison, we also plot the widely used analytical fit of [Sheth] 


et al. (2001) at each redshift which has been tested against 


a wide range of numerical simulations. The halo mass func¬ 
tion of our different simulations are in good agreement with 
the analytical model over the mass range where agreement 
should be expected. The minimum mass of haloes assumed 
to host ionising sources is indicated by the vertical dashed 
lines as well as the corresponding circular velocities. 


APPENDIX B: THE LUMINOSITY FUNCTION 
OF IONISING SOURCES 

In our ATON simulation we have assumed that the lumi¬ 
nosity of ionising sources scales linearly with the mass of 
the dark matter haloes identified in the RAMSES simu¬ 
lations as discussed in Appendix [A] As discussed by e.g. 


Trenti et al. (2010) models with a linear relation between 


dark matter halo mass function and galaxy luminosity fit 
the evolution of the high-redshift luminosity function rea¬ 
sonably well. To test for consistency of our modelling and re¬ 
late our assumed population of ionising sources to observed 
high-redshift galaxies we show here the evolution of the lu¬ 
minosity function obtained from such a linear relation be¬ 
tween galaxy luminosity and halo mass and the halo mass 
functions discussed in Appendix [A] 

In our source model, the total comoving ionising emis- 
sivity is scaled to that i n the HM2012 UV background model 
as described in section 2j[ From this, we can calculate the 
integrated emissivity £ 1500 (z) at 1500 A at a given redshift. 
Summing over all ionising sources gives then the normaliza¬ 
tion between halo mass and 1500 A flux /isoo- 

Figure [Bl] compares the resulting luminosity functions 
with the observed luminosity function of high-redshift galax¬ 
ies. At the highest redshifts with available data z ~ 9 — 8, 
there is excellent agreement of our simple model of the lu¬ 
minosity function and the data. At lower redshifts we some¬ 
what overpredict the faint end of the luminosity function. 
Note, however, that we are interested in the hydrogen ionis¬ 
ing luminosity for our reionisation simulations and that the 
escape fraction of ionising photons is generally believed to 
increase with decreasing galaxy luminosity. The somewhat 
steeper luminosity function for the ionising flux we have as¬ 
sumed is therefore probably as good a bet as currently pos¬ 
sibly. 


APPENDIX C: MEASURING THE MEAN 
FREE PATH OF IONISING PHOTONS FROM 
SIMULATIONS 


The mean free path of ionising photons in simulations can be 
computed from the remaining transmission after a specific 
path length and by assuming that the average transmission 
decreases exponentially with increasing distance. A more so¬ 
phisticated method is to adjust the length scale over which 
the transmission is measured depending on the value of the 
mean free path (Emberson et al. 2013). 

Here we use a more direct approach which does not rely 
on choosing a length scale for the measurement or on assum¬ 
ing a functional form of how the transmission decreases with 
distance. Instead we directly use the definition of the mean 
free path, 


Amfp — j d/ >• 


(Cl) 


Here f(x) = exp(— r(x)) is the transmitted ionising flux at 
path length x and (...) denotes an average over many lines 
of sight. The only assumption we need is that the neutral 
hydrogen density is constant within each simulation cell. 

Without loss of generality we can assume that the initial 
ionising flux f(x = 0) = 1. In this case the integral in the 
denominator in Eq. ( |C1| is equal to —1, so that we are left 


1 The comoving emissivity table of HM2012 
(in erg s -1 Hz -1 Mpc -3 ) can be downloaded at 

http: //www.ucolick. org/~pmadau/CUBA. 
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Figure Al. Evolution of the halo mass function in our different simulations. The analytical model of[Sheth et al. (2001) is shown by the 
black solid curves. The coloured vertical lines show the minimum mass assumed to host ionising sources in each model. The corresponding 
halo minimum circular velocities are displayed on the plot (in km/s). 


with 


X mfp = -{[° xdf) = (JfX i ), 


(C2) 


where A i is the contribution of the z-th cell in the line of 
sight to the integral. Next, we note that 

- df = e~ T(x) ^ dr = e~ T(x) ff- dx. (C3) 

J Ax Ax v ' 

In the second equality, we have used the assumption that 


the neutral hydrogen density is constant within each cell so 
that dr/dx = n/Ax, were r* is the optical depth of cell 
z and Ax is the cell size. Before performing the integration 
over cell z, we change the integration variable by defining l = 
x — (i — 1) Ax, which is just the distance from the beginning 
of cell i. Using this we find 


r(x) = J2rj + 

3 = 1 


(C4) 
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Figure Bl. Evolution of the luminosity function in our different simulations assuming a linear relation between halo mass and 1500 A 
flux as described in the text. Black and purple squares with errorbars show respectively observational constr aints fromfF inkelstein et al.| 
(|2014|> and|Bouwens et aT] ( |2014| >. The solid curves are fits to observational data from |Bouwens et al.| ( |2007| ), |Bouwens et al.| ( |2008| > and 
|Trenti et al.| ( |2010| ). 


so that we can write A^ as 

_ /* A * 
Xi = e ^o= 


lTj / [(i - l)Ax + l\ e A 

Jo 

— e ~ ^i=i Tj [(1 — e~ Ti )(i — l)Ax 
1 - e~ Ti (1 + n) 


- 

- ■o = 

1 


l 

Ax 


+ Ax- 


Ti 


(C5) 


Using this expression, we find the mean free path of ionising 
photons along a specific line of sight by summinng over the 


cells along it, i.e. by computing JA We perform this sum¬ 
mation at least until we have reached the opposite side of 
the simulation box. If the transmitted flux fraction there is 
larger than 10 -8 , we continue the summation along another 
randomly chosen line of sight. We go on in this manner until 
f(x) < 10 -8 . Finally, we repeat this prodecure many times 
for different lines of sight and average the obtained mean 
free path values to get a global estimate of the mean free 
path of ionising photons. 
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Figure Dl. Right: The photoionisation rate Ti 2 due to bright sources at redshift z ~ 5.8 for a slice of our bright source model assuming 
a mean free path of A^ f p ~ 35 comoving Mpc/h and 2 bright ionising sources in the periodically replicated (100Mpc/h) 3 volume. The 
photo ionisation rates are normalised to match the observed cumulative Lya optical depth as described in the text with /3 ga i = 1.15 
and Alight = 1-^5 and the mean photoionisation rate indicated on the plot is averaged over the whole simulation box. Left: The 
photoionisation rate Ti 2 in the ‘512-20-goodT ATON radiative transfer model, also at redshift z ~ 5.8. 


APPENDIX D: A TOY MODEL FOR RARE 
BRIGHT GALAXIES/(FAINT) AGN 


In section |3.3.4| we have argued that fluctuations of the 
photoionisation rate caused by rare bright galaxies or per¬ 
haps more likely (faint) AGN may be responsible of the 
broad Lya opacity PDF measured for 50 Mpc/h chunks at 
z ~ 5.6 — 5.8. Here we present details of the toy model dis¬ 
cussed there that we have used to test the effect of rare 
bright sources on the Lya opacity PDF at redshift z ~ 5.8. 

We use the location of the most massive dark matter 
haloes in a periodically replicated large volume simulation 
with a box size 100 Mpc/h on a side (1200 3 particles, also 
used and described in Choudhury et al. 2014) at z ~ 5.8 as 
positions for the bright ionising sources. We then compute 
the photoionisation rate T due to these bright sources at 
every position in the 100 Mpc/h volume with the simple 
attenuation model used by Becker et al. (|2015| (their ga laxy 
UVB model of section 4.2, but see also |Bolton &; Vie1||201 1| 
and Viel et ah]|2013 ). 

We assume a spectral energy distribution appropriate 
for AGN for the bright sources of the form (see |Vanden 
Berk et al.|2001 and Telfer et al.||2002 ), 


L v oc 


z/- 0 44 if A > 1300 A 
V - 1 - 57 if A < 1300 A 


(Dl) 


As a first guess we then calculate the L u ( 912) luminosity for 
each source such that the total eg 12 comoving emissivity in 
the v olume is equal to the one derived by Giallo ngo et al.| 
(20 1 5), 2.5 x 10 24 erg/s/Hz/Mpc 3 . This results in a total ion¬ 
ising emissivity of 5.18 x 10 56 and 3.63 x 10 56 photons per 
seconds for the two sources. 

At each spatial position, we compute the specific in¬ 
tensity of the ionising background between 1 and 4 Ryd by 


summing over the contribution from each bright source, 

N -j- , \ -|r,;-r| ( u ^-3(/3-l) 

Li{n,u) A^ f | U912J 




“ 47r|r* - ] 

Z=1 1 
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where, V 912 is the frequency at the HI ionising edge, and 
/3 = 1.3 is the slope of the HI column density distribution, 
which gives the dependence of mean free path on frequency 
(see also Son gaila &; Cowie|2010 and Becker &; Bolton|2013 ). 
The sum in equation |D 2 1 is performed for all sources within 
the periodically replicated simulation boxes (in practice, we 
included only the sources within the 100 Mpc/h box and its 
26 directly neighboring periodic replications). We assume a 
mean free path A^f 2 as extracted from our ATON ‘512-20- 
good_h’ full radiative transfer simulation which has a value 
35.09 comoving Mpc/h at redshift z ~ 5.8. 

The HI photoionisation rate due to bright sources is 
then computed as 


of AS& 
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r(r) = 4 tt / 
J Vi 


41/912 dv 

^912 ' lV 


(D3) 


where (Jui{y) is the photoionisation cross-section (cal¬ 
culated from the fit of jHui & Gnedin 1997). 

We combine the photoionisation rates due to bright 
sources in our toy model as follows with the ionising UV 
background due to the much more numerous galaxies driv¬ 
ing reionisation in our ATON simulations to calculate the 
expected effect of bright sources on the Lya opacity PDF. 

• We randomly choose a line-of sight through the 
100(Mpc/h) 3 volume (along one of the principal axis) for 
which we have modelled the contribution of bright sources 
to the photoionisation rate Tbright- 

• We concatenate three randomly selected line-of-sights 
from the ‘512-20-good_l’ full radiative transfer simulation 
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and place them along the line-of-sight through the bright 
source model and call the photoionisation rates due to the 
sources in the radiative transfer simulation r ga i. 

• We calculate the combined photoionisation rates along 
the line-of sight, Fg a l _|_ bright — /^galTgal ~\~ /^brightTbright • 

• The value of /bright is calculated such that for chosen 
values of the number of bright sources and /3 ga \ (the two ad¬ 
ditional free parameter of our model once a mean free path 
has been chosen), the mean value of the total photoionisa¬ 
tion is that required to best match the Lyct opacity PDF. 
Note that when we calculate (r ga i + bright), we do so for the 
entire 100 Mpc/h cube in order to have the same value of 
/^bright in each line of sight. 


For each model for the spatial fluctuations of the pho¬ 
toionisation rates we calculate 5000 mock Lyct spectra for 


50 Mpc/h chunks as in the observed sample of Becker et al. 
( |2015 ) using the density, temperature, and peculiar velocity 
fields from the ‘512-20-goodT model. We have then var¬ 
ied the number of bright sources and /3 ga \ with the aim of 
matching the observed Lya opacity PDF. In practice for a 
chosen number of bright sources we have started with choos¬ 
ing /3 ga i = 1 and determined the value of /bright required to 
match the mean photoionisation rate in the ‘512-20-good_h’ 
model at z = 5.8 which is reasonable close to that required 
to mach the observed Lya opacity PDF. We have then fur¬ 
ther moderately rescaled r ga i + bright to match the observed 
Lya opacity PDF as well as possible. 

The model discussed in section 13.3.41 for which we were 
able to obtain a good match had two bright sources, /3 ga \ = 
1.15 and /bright — 1-25, but there will certainly be other 
parameter choices that will equally well reproduce the data. 
In this model the galaxies driving reionisation and the bright 
sources responsible for the large scale fluctuations contribute 
about equally to the integrated ionising UV background. In 
figure [ET] we show the photoionisation rate Tbright in a slice 
of the 100(Mpc/h) 3 volume for this model. For comparison, 
we also show a map of a slice of r ga i. As expected our toy 
model produces of order unity fluctuations of T on scales of 
> 50 Mpc/h, while the spatial fluctuations due to the much 
more abundant galaxies in our radiative transfer simulations 
occur on much smaller scales. 

It should, however, be kept in mind that the model pre¬ 
sented here neglects many aspects which will be important 
for an accurate modelling of the effect of bright sources. 
We have e.g. not self-consistently modelled the effect of the 
bright sources on the mean-free path of ionising radiation. 
The large spatial extent over which individual bright sources 
dominate also means that light-travel time effects become 
important if the sources are short-lived as is likely for bright 
starbursts and AGN. The angular distribution of the emit¬ 
ted ionising radiation is also likely to be not isotropic. We 
should thus emphasise here that the results regarding the 
Lya opacity PDF based on this simple model should just be 
seen as proof of concept. 







